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1.0 Introduction 



The subject of linear optimal approximation has received considerable 
attention in recent years [4], [5], [6], [7], [8]. The subject of multi- 
variate approximation for scattered data, including optimal approximations, 
is reviewed in [9]. The idea is appealing since the optimal approximation 
in a certain space of functions minimizes the norm of the error functional 
for approximations in that space. When the space is a Hilbert space, the 
computation of optimal approximations becomes rather simple, in theory [2], 

A known reproducing kernel function provides the representers of linear 
functionals defined on the space. The optimal approximation satisfies the 
system of equations obtained by requiring that the approximation be exact 
for the representers of the functionals being used for the approximation, 
usually point evaluation functionals. 

In practice, it seems that optimal approximations have not been used 
very much. This is perhaps partly because of a lack of experience with them, 
as well as the fact that use of the representers as a basis set for optimal 
approximation is the analog of the use of truncated power functions as a 
basis set for univariate spline approximation. 

The particular space of functions of two variables to be considered here 
are the Sard "corner spaces". Bp -j [8]. A suitable completion of these 
spaces into a Hilbert space and construction of the reproducing kernel was 
recently accomplished [1]. Since these spaces are made up of functions whose 
partial derivatives, up to a certain order in each variable, are absolutely 
continuous, these spaces contain spline functions in two variables, and the 
optimal approximations are splines. The Sard corner spaces have the property 
that the representers reduce to products of functions of one variable, thus 
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a base point (a,b) in the region of interest. While this point can theo- 
retically be anywhere, practically it is desirable for it to be at one of 
the corners of the region of interest (assumed rectangular). There are two 
reasons for this; (i) Computation of the representers is simplified; and 
(ii) the representers have continuous partial derivatives in x only up 
through order p - 1 at x = a and in y only up through order q - 1 at 
y = b , while at other points partial derivatives in x and y are contin- 
uous up through order 2p - 2 and 2q - 2 , respectively. The second reason 
is the primary reason we assume that (a,b) is the origin and that we are 
only interested in (x,y) points in the first quadrant. 

In connection with the previous paragraph, we note that approximations 
in Bpp -j are not invariant with respect to translation (unless the base 
point (a,b) is also translated), nor with respect to stretching or shrink- 
ing of the coordinate system. The author had previously commented that they 
were invariant [3]. The base point (a,b) is the point at which Taylor 
series (with remainder) for the functions exist, and the approximations are 
clearly dependent on that point. 

In the general case, the reproducing kernel function for -j is of 

the form K(a,b;u,v,x,y) = gp(a;u,x)gq(b;v,y) where for a < u,x , 




i<p 




Here the notation w^ 



means w / i ! while 



w , w > 0 



w 



'+ 



0 , w < 0 
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For the case we consider, with a=b=0, p=q , and where the func- 
tionals are point evaluations, say at the point (x^.y^) , the representers 
are K. (x,y) = g (0;x. ,x)g (0;y . ,y) . As in what follows we have 

i H H 

here used simplified notation for the representer associated with (x..,y.) . 
We will consider q = 1 and q = 2 and thus note that g^ (0;x^ ;x) = 1 + x - 
(x - x i ) + and g 2 (0;x i ;x) = 1 + x^x + ^ x^x 2 - ^ x 3 + ^(x - x..)^ . 

One observation about optimal approximations is now in order, and should 
lead to increased interest in their use for approximations where the data is 
irregularly spaced. The data points themselves generate the representers, 
and hence a set of basis functions for the approximations. One does not 
need to be concerned about whether or not the basis functions have the inter- 
polation property on the set. Unlike more common basis functions, e.g., 
polynomials, the representers naturally form a linearly independent set over 
the data points. This is not meant to imply that the representers of point 
evaluation functionals are well suited to computation, however. We treat 
this problem in more detail in later sections. 

1 .1 The interpolation problem 

The underlying problem we shall be considering is that of function 
approximation by interpolation for functions of two or more variables. The 
case of more than two variables is a straightforward, if tedious, general- 
ization, and the discussion is limited to two independent variables. Assume 
that the points (x^.y^jZ^) , k = 1,...,N are given. No assumptions are 
generally made as to the spacing of the points, although in some instances 
we will consider special cases. We assume that if i * k, (x^ ,y.) * (x^y^) 

Approximation problems other than interpolation are treated in identical 
fashion. One obtains the same coefficient matrix for approximate integration 
or differentiation, for example. 
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2.0 Optimal approximation in -j, 

In this section we will consider in some detail the problem of computing 
optimal approximations in the Sard space B ^ ^ . Again we emphasize that 
the region of interest is assumed to lie in the first quadrant and that the 
base point (a,b) is taken to be the origin. 

2. 1 Representers of point evaluation functionals as a basis 

As noted in the introduction, the representers of point evaluation 

functionals have the form K.(x,y) = [1 + x - (x - x.) + ][l + y - (y - y.) ] . 

J J J 

These functions K. are continuous, with first partial derivatives in x 

J 

and y which have jumps at x = x. and y = y . , respectively. In each 

J J 

of the rectangles [0,x .]x[0,yj] , [0,x ,.]x[y . ,») , [x^ ,»)x[0,y .] , and 
[x • ,°°)x[y . ,=°) the function is bilinear. In the latter rectangle it takes 

J J 

on the constant value (1 + x.)(l + y.) . 

J J 

Use of the K.(x,y) as basis functions appears to be suspect. The 

J 

coefficient matrix (the Gram matrix) is symmetric, but casual observation 

would lead one to suspect it is not particularly well conditioned. We shall 

see that it is better than the author's inclination toward it. The system 

of equations has the form 
N 

(1) 2 A .K . (x . ,y . ) = z,, i = 1 , . . • ,N . 

j=! J J 1 1 1 

Some numerical experiments were conducted to compute the condition number 
(with respect to max row sum norm) of some Gram matrices of various sizes. 

Points were generated by a random number generator in the square [0,10]x[0,10] , 
the Gram matrix was formed, and the condition number computed. In case (i) the 
points were allowed to be anywhere in the square. In case (ii) the square was 
subdivided into [/N] squares, then one point was generated at random in each 
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smaller square, and any remaining points were generated at random throughout 
[0,1 0]x[0,l 0] . In cases (iii) and (iv) the points in (i) and (ii) were 
translated to the square [90,1 00]x[90,l 00] . In cases (v) and (vi) the 
points in (i) and (ii) we stretched over the interval [0,1 00]x[0,l 00] . 

In case (vii) points with integer coordinates in [0,10]x[0,10] were selectee 
by a random number generator. These cases give some indication of condition 
numbers and the effects of translation and stretching. No claim is made that 
the cited results are representative nor that an exhaustive set of calculatior 
were made. Other computations, not tabulated here, follow the same trend, 
however. The results of the calculations are tabulated in Table 1. Mot all 



combinations were 
Case N -> 


computed. 

10 


25 


50 


100 


(i) 


170 


3400 


6700 


67000 


(ii) 


420 


1500 


5700 


36000 


(iii) 


24U0 


48U0U 


18UUUU 




(iv) 


3600 


37000 


200000 




(v) 


190 


3200 


8600 




(vi ) 


600 


6100 


6200 




(vii) 




2900 


15000 





Table 1: Condition Numbers of (Kj(x^,y^)) for Bp ^ 

One can make several observations from the table. First, for moderate 
values of N , even up to 50 or 100, satisfactory results can easily be obtain 
by computing with the representers as basis functions. Depending on the 
accuracy required in the computed answer and the precision one can (or is 
willing) to use in the computations, N could be quite a large number, per- 
haps as large as is feasible to consider for a global approximation. Second, 
the effect. of stretching on the condition number is rather mild. Third, the 
translation of the points away from the base point increases the condition 
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number, indicating that one should select the base point (a,b) close to a 
corner of the region of interest, as we have done. 



2.2 A set of basis functions with compact support 

It is possible to construct a set of basis functions with compact support 
which leads to a block tridiagonal system of equations, where each block is 
tridiagonal, for approximations in Bj-^ ^ . Unfortunately the scheme is 
practical only in the special instance that the data points lie on a grid. 

Such a grid exists for any set of points, but ordinarily only one point lies 
on each horizontal and each vertical grid line. For the scheme to be prac- 
tical it is not necessary that each grid point be a data point, but many 
points should be on most grid lines. This will be made more explicit after 
the development. 

We must alter our usual notation slightly. Suppose we have grid points 
(x. ,y-)» i = 1 , n, j = 1 m . We assume that the x. and y. are in 

I J * J 

increasing order. Denote the set of subscript pairs corresponding to data 

points by I . Then corresponding to each (k,z)el there is a known function 

value z k . Nondata points on the grid will be denoted as ( i , j ) / I , where 

we will always assume that 1 < i < n and 1 < j < m . 

It is easy to obtain functions of one variable of the appropriate form 

which have compact support. They are 

n 

G • (x) = E a. g n (0;x ,x) , i = l,2,...,n , 
i r=1 ir r 

where as before g^ (0;x r ,x) = 1 + x - (x - x^) + , 

1 + x„ 



with a 



11 = (1 + x,)(x, - xj ’ a l 2 



1 



a . 



^ -j f \ A ^ XX 1 j It- A -j A ^ 

1 X . , . - x. -i 

1+1 1-1 



ii-1 



x i-l ‘ x i 



• a ii ■ (x, - x 1+1 Hx 1 _ 1 - x,) * “11+1 ' x, - x 1+1 



- 9 - 



1 



1 



""- 1 Vl - x n 



nn x - x , 
n n-1 



and a. =0 if |i - r| > 1 . For convenience, let x n = 0 and x 

ir 11 0 n+1 



then we note that G. (x.) = 6.. as well as G. (x) * 0 only over the interva 

* J ' J * 

(x. i , x^ + -|) • Also construct the dual functions Hj(y) , with 

Hj(y) = 2 3j S g-|(0;y s ,y) . Note that the G i and HL are linear B-splines. 



To satisfy the interpolation requirements in terms of the local basis 

n m 



1 1 III 

functions G. (x)H .(y) ,we obtain the equations Z z a. .G.(x. )H .(y ) = 
1 J i=l j=l 1J 1 K J * 



a, = z 



k,£ k,£ 



, (k,A) el. For (k,£) £ I , the products g ] (0;x k ,x)g 1 (0;y £ ,y 



cannot appear in the approximation. Substituting for G.(x) and H.(y) the 

3 



n m 



1 1 m 

approximation becomes Z z a-. G.(x)H.(y) = Z a., z a- g n (0;x ,x) z 

i=l j=l 1J 1 J i,j 1J r lr 1 r s 



0 js 9^ 0 »y s ’ y) = - Z . a i j Z a ir 0 js g l (°’ x r »x)g 1 (0;y s ,y) 

t >J r S s 



We 



then set the coefficient of g-j (0;x k ,x)g-| (0;y, ,y) for ( k , Jt) i I 



equal to zero, obtaining as the system of mn equations for the a.. , 

1 3 



( 2 ) 



f 



a M “ Z M ’ (M) e 1 



) 

V 1,J 



a ij “ik 6 j* = °> < k - z > *■ 1 ■ 



If we order the equations and variables in some logical fashion, say 
(1,1),.., (l,m), (2,1),..., (2,m),..., (n,m), for m < n , the fact that 



‘ik 



= 0 if | i - k | > 1 and 0. = 0 if |j - a| > 1 shows that the 

3 ** 
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resulting coefficient matrix of the system of equations is of the form 



T 



11 

r 21 



\ 



'12 

r 22 



23 



i •> 

nn- I 



0 \ 



nn 



Where each T. . = a. . 

"1 J "1 J 



e 



n 



e 12 o 



>21 



22 



23 



0 

0 



mm-1 



mm 



/ 



except that a row of T corresponding to (k,«,) e I is replaced by a new 
row with zeros everywhere except for the unit diagonal element. 

We can note that ar| d ^jj+1 are non P osl ti ve and except for 

j • 1 and j ■ m. 6 jd • - Sjj., - B jJ+1 , while for j = 1. 6,, > - e, 2 

and for j = m , = - 3„ m •, . Thus each block T.. is diagonally 

mm mm-1 ij 3 J 

dominant. An interesting aside is that the system (2) looks very much like 
the system of equations obtained when solving Laplace's equation by finite 
differenceson a rectangular grid, where a 9 point approximation to the Laplacian 
is used. Here, of course, the "boundary conditions" may be scattered through- 
out the region. In this respect the coefficient matrix is well suited to the use 
of iterative methods for solution of (2). 

2 

For the usual case, the above leads to a system of N equations, and the 
amount of computation will be too large for even moderate N. Use of a block 
elimination equation solver reduces the problem to one of repeated solution 
of n systems of equations of m equations (usually full, but not always). 



- 11 - 



3 

Although the original problem can be solved in about N /6 operations, the 

new basis requires about n m^ , or about operations if m = n = N . 

One should point out, however, that after obtaining the a., one has sufficie 

information to obtain the value of the interpolating function by bilinear inte 

polation since a^ represents the function value at (x^ ,yj) . 

We will consider the case in which there are several data points on each 

grid line. The amount of work required to solve (2) decreases somewhat as 

the fraction of known grid values increases, however, the major effect is that 

the relationship between (n,m) and N changes. Assume that the fraction 

of grid points at which data is known is p, 0 < p < 1 . The total number 

of data points is then N = pmn , and the number of operations required to 

3 3 3 

solve (1) is about (pmn) /6 . For (pmn) /6 a nm we see that 

q o 3 o vq 

p J ~6/n^ , or p a x\~ c . Representative values of these fractions are 
given in Table 2, along with the total number of equations in (1) and (2). 

If p is larger than the listed value in Table 2, (2) can be solved in fewer 
operations than (1 ) . 



n 


P 


N 


mn 


4 


.721 


2.88m 


4m 


8 


.454 


3.63m 


8m 


15 


.299 


4.48m 


15m 


50 


.134 


6.69m 


50m 


100 


.0843 


8.43m 


100m 


1000 


.0182 


18.2m 


1000m 



Table 2: Fraction p of grid points to be 
known for comparable number of operations in 
solving (1) and (2). 
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2.3 Basis functions which are zero in some regions 



Because the representers of point evaluation functionals are constant for 
large enough values of the independent variables, a certain linear combination 
of any two can be made zero for large enough values of the independent vari- 
ables. In particular, consider K.(x,y) and K.(x,y) . The function 

' J 

(1 + x. )(1 + y.j)Kj(x,y) - (1 + Xj)(l + yj)K..(x,y) is zero for x > max^. ,Xj) 

and y > max(y.. ,yj) . Thus a new set of basis functions with zero values over 

part of the region of interest can easily be constructed. It is desirable to 

first order the data points in terms of their "distance" from the origin. It 

seems reasonable to order the (x^^y^) in terms of nondecreasing values of 

(1 + x k )(l + y k ) . This is not the only ordering which can be used, but it 

t h 

carries the assurance that the new j basis function will be nonzero at the 

x. L. 

j point, and also has an added benefit we will discuss later. 

Assume that the data points are ordered so that p = (1 + x )(1 + j .), 

K K K 

k = 1,...,N is a nondecreasinq sequence. 



Then define L.(x,y) 
J 



- 0jt ,K.( x y) 

D /j + i (x j' y j) - VlVVV 



and L N (x,y) 



K N (x,y) 

V X N’V 



j = 1 



,N 



1 



Then the Lj(x,y) satisfy L j ( x ^ ,y.) = 1, j = 1,...,N and 
Lj(x,y) = 0 for x > x^ = max(Xj ,x^ +1 ) and y > y. = max(yj ,y^ +1 ) , j = 1,...,N- 

Using these basis functions, the interpolation problem requires the so- 
lution of the system 
N 

(3) l A.L,(x.,y.) = z., i = 1 N . 

j=l J J 1 1 1 



- 13 - 



Because of the above ordering on the data points, there is some possibility 



of the entry L-(x.,y.) being zero for j < i , that is an element below 

J * 



the diagonal of the coefficient matrix. 



The basis functions L.(x,y) have a property which is desirable and 

J 



which arises out of our ordering of the data points. 

Proposition : In the first quadrant, |L.(x,y)| < 1 . 

J 

We note that it is also true that if Pj + -j < P j , then |L-(x,y)| < 1 



except along one of the rays which start at (x.,y.) and extend horizontally 

J J 



to the right and vertically upward. The proof is simple and will not be 



given. It consists of considering the value of L.(x,y) at all points where 

J 



the first partial derivatives are discontinuous, since any extrema must occur 
at such a point. These points are (0,0), (xj ,0), (x^.O), ( 0 ,y ^ ) , (0,y. +1 ), 
(Xj ,y .), ( x j ,x j+l ) » ( x j+i>yj) » and (x- +1 ,y . +1 ) . Typical function values 



are shown in figure 1. Recall that in each rectangle L.(x,y) is bilinear an 

J 



hence determined by its values at the corners. 
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-.286 


-.857 


0 


-.286 


-.857 


0 




( x 2 » y 2 ) = ( 2 > 8 ) 




.048 


.143 ^ X 1 , y l ^ = ^ 6s2 ^ 


1 .000 


.016 


.048 


.333 



0 2 6 10 



Figure 2: Values of L-|(x,y) 



0 



0 



.000 



.333 
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Certain behavior can be classified further. For example: (i) if 

Xj + 1 > x j ,y J+1 > Yj , then |Lj(x,y)| < 1 except at (x^.y.) ; (ii) 

if p j + i = p j , then Lj(x,y) = 0 for x <_ min(x^ ,Xj +1 ) and 

y < min(y.,y. +1 ) and L j( x j+i »y J+ i ) = -1 . 

The coefficient matrix of the system (3) has its largest element (in 

magnitude) on the diagonal, and some zeros may occur below the diagonal. 

The data points used to generate Table 1 were used to test the effect- 
iveness of the introduction of zeros and to determine the condition numbers 
associated with the new basis functions. The results are shown in Table 3, 
with the number in parenthesis indicating the number of 1 eadinq zeros in 
the matrix. By reordering the columns of the matrix it is sometimes possible 
to introduce many more leading zeros, and while a scheme of this sort has not 
been implemented, in many cases it would substantially reduce the number of 
operations required for solution of the system (3). 

The system (3) does not have a symmetric coefficient matrix, and unless 
approximately 30% leading zeros are introduced, it will take fewer operations 
to solve (1) than to solve (3). However, the condition number of the new 
coefficient matrix has been smaller in every case examined, sometimes by a 
factor of 15 or more, but more commonly by a factor of 2-5. Reordering 
columns for a maximum number of leading zeros will often result in 30% or 
more leading zeros, based on some hand computations where no effort was made 
to obtain the maximum number of leading zeros. 



- 16 - 



Case N -* 


10 


25 


50 


100 


(i) 


36(4) 


250(20) 


1800(261) 


18000(535) 


( 11 ) 


32(18) 


460(32) 


2000(162) 


12000(582) 


( i i 1 ) 


1800(6) 


2600(39) 


30000(315) 




(iv) 


160(18) 


3300(62) 


19000(288) 




(v) 


23(2) 


150(20) 


1400(195) 




(vi) 


28(18) 


390(16) 


1500(105) 




(vii) 




620(21) 


6200(146) 





Table 3: Condition Numbers of (L^ (x - ) ) and Number of 

Leading Zeros (in Parenthesis) 



3.0 Optimal Approximation in Br 2 2 ~| • 

The problem of computing optimal approximations in the Sard space 
B[2 2~| 1S somew ^ a t more difficult than for ^-j . The reproducing kernel 
functions are seen to be piecewise bicubic functions, reducing to bilinear 
functions for sufficiently large values of the independent variables. We 
shall investigate the feasibility of extending the results of the previous 
section to Br 2 2 -j this sec tion. 



3.1 Representers of point evaluation functionals as a basis 

The representer of the point evaluation functional at the point (x.,y.) is 

2 3 ^ 3 

..(X,y) = [l + y + 2 Xj X - f X 3 + 5 <X - X l)+ }[1 + y.y + | y. 2 y - ^ y 3 X - y.) + 

This function is cubic in x for 0 < x < x. , and linear in x for 

J 

x > Xj , and the dual holds in y . These functions increase rapidly since the 

point x = x. is the inflection point of the cubic in x , and thus when the 
3 
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function is linear in x , it has slope the same as the maximum slope of the 
cubic. Because of this, the Gram matrix is not well conditioned. Following 
a similar path to that taken in the preyious section, some condition numbers 
for the Gram matrix were computed for some sets of randomly generated points. 
The results are given in Table 4. The point description column refers to the 



descriptions in 


section 2.1 


• 






Case N ■* 


10 


25 


50 


100 


(i) 


3.94-10 4 


3.5M0 7 


8.27-1 0 7 


4.58-10 


(ID 


3 . 89 - 1 0 5 


1.89-10 6 


4 . 47 - 1 0 9 


9.95-10' 




Table 4: 


Condition Numbers of 


(K j (x i ,y.) for 





The observations which we wish to make are that: (1) One will quickly be 
in numerical trouble in Real *4 on the IBM 360, and (2) While the condition 
number is large, meaningful computations can be done in Real *8. 



3.2 A set of basis functions with compact support 

A similar construction for Bp, 2 "j as was pursued for Bp p in section 

2.2 can be done. There is some question as to what conditions should be imposed 

on G-j , G 2 , G 3 , G 4 , G n _ 3 , G n _ 2 , G^_-| , and G^ , (as well as the corresponding 

dual functions H. ), but several reasonable options are open. In the general 

n 

case one wants G^x) = £ a. g 2 (0;x ,x) so that G- (x..) = 1 , G^(x) = 0 , 



x <_ x.j _2 or x > x - + 2 • Proceeding in similar fashion one will obtain a system 
of equations which in block form has 5 non-zero blocks, each block being a 
square matrix with 5 non-zero elements per row. Numerically this is somewhat 
more complicated than before, but there are instances where it could be useful. 
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No details of the construction have been carried out here, but the general- 
ization is straightforward. Also note that the G, and H. here are cubic 

* J 

B-spl ines. 

3 . 3 Basis functions which are zero in some regions 

It seems natural to be able to extend this idea to B^ £-| • The basis 

functions K.(x,y) are bilinear for x •> x . and y > y. , thus it seems 
J J J 

possible a certain linear combination of them could be made identically zero 

to the right and above all points (x.,y.) associated with those five basis 

J J 

functions. Proceeding in the obvious fashion, ordering the points ( ,y ^ ) 
by some rule, we then wish, for j <_ N - 4 , to construct functions 
j+4 

Lj(x,y) = Yj k K k (x,y) such that 
k-J 

Ljfx.y) 0 for x > x% max (x j+4> x j+3 , x J+2 , x j+ , , Xj) 

and y > yj * = max (y J+4 . y j+3 , y J+2 , y J+1 , yj ) 

and Lj(xj,y .) = 1 . Unfortunately, if the points (x j+4 ,y j+4 ), (x j+3 ,y j+3 ) , 
( x j +2 >yj +2 )> ( x j + i »yj +] ) li e on any bilinear curve, this system of equations 
generally has no solution. Thus the ordering imposed earlier would at least 
have to restrict one away from four successive points lying on a bilinear 
curve. In general, this is not possible. 

Because of the rather more restrictive region where the function is zero, 
less benefit is likely to accrue anyway. In addition, even when one can con- 
struct such sets of basis functions, it is not possible to bound the function 
Lj(x,y) as in the Bp ^ case, and in particular it cannot be bounded by one. 
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4.0 Conclusions 



This investigation has determined that use of the representers of point 
evaluation functionals can be used as a basis for some problems without en- 
countering severe computational problems. For smooth approximations this is 
probably not generally true, however. In addition, because the computational 
burden for global approximations is likely to be quite large, it is the 
author's opinion that local approximations must be investigated for smooth 
interpolation. The time is perhaps propitious for an investigation into the 
underlying mathematical basis for some previously suggested schemes for local 
smooth interpolation. 

The use of optimal approximations in Sard corner spaces for the inter- 
polation of irregularly spaced data results in an approximation which has 
discontinuities along the lines parallel to the axes through each data point. 
This would seem to the author to be unnecessarily complicated, and the author 
intends to investigate global approximations in which the discontinuities 
are less numerous and local in character rather than extending alona lines to 
infinity. 
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